Complex SIR model with Variability

Complex SIR Model

Assumption

• The recovered (R) compartment will no longer be susceptible

Parameters with default values:

• Population \((N = 300)\)

• Transmission/Infection Rate \((r = 0.02)\)

• Percentage of healthy individuals with symptoms \((g = 0.2)\)

• Initial number of infected individuals day1 \((d_1 = 5)\)

• Initial number of cured individuals day1 \((c_1 = 0)\)

• Percentage of individuals receiving Treatment A \((PTA = 1)\)

• Percentage of individuals receiving Treatment B \((PTB = 0)\)

• Treatment A effectiveness \((A_{eff} = 0.75)\)

• Treatment B effectiveness \((B_{eff} = 0.25)\)

• Treatment A cost \((A_{cost} = 20.00)\)

• Treatment B cost \((B_{cost} = 15.00)\)

• Number of days since first case \((day = 38)\)

Model upgrade with randomness - Additional Variables

• Initial number of known cases (untreated) day1 \((u_1 = 0)\)

• Percentage of individual being tested daily \((PS = 1)\)

• Test sensitivity \((tp = 0.99)\)

• Test specificity \((tn = 0.95)\)

• Testing cost \((Test.cost = 1)\)

complex.model.rand = function(N=300, r=0.02, g=0.2, d1=5, c1=0, u1=0, PS=1, tp=0.99, tn=0.95, PTA=1, PTB=0, A.eff=0.75, B.eff=0.25, day = -1, A.cost = 20, B.cost = 15, Test.cost = 1){
  
  # Initial conditions
  curr.day = 1
  d = d1
  c = c1
  u = u1
  # Calculate initial statistics
  h = N - d - c
  hs = round(g * h)
  s = hs + d
  a = s - u
  ns = round(PS * a)
  ad = round(PS * (d - u))
  as = ns - ad
  tps = rbinom(1, as, 1-tn)
  tpd = rbinom(1, ad, tp)
  tpt = tps + tpd
  tpt1 = tpt
  NDAT = tpd + u
  NDTA = 0
  NDTB = 0
  DCA = 0
  DCB = 0
  NAT = 0
  NTA = 0
  NTB = 0
  SDC = 0
  New = min(N - c - d, round(r * hs * d))
  
  # Table storing S, I, R compartments over time
  df = data_frame(Day = curr.day, 
                  Healthy = h,
                  Symptopms = hs,
                  Diseased = d,
                  Tot.Symp = s,
                  Avail.screen = a,
                  Untreated = u,
                  Screen.Tot = ns,
                  Screen.Sym = as,
                  Screen.Dis = ad,
                  Test.Pos.Sym = tps,
                  Test.Pos.Dis = tpd,
                  Test.Pos.Total = tpt,
                  Avail.TRT.Dis = NDAT,
                  TRT.A.Dis = NDTA, 
                  TRT.B.Dis = NDTB, 
                  TRT.A.Cure = DCA, 
                  TRT.B.Cure = DCB,
                  Avail.TRT.Tot = NAT,
                  TRT.A.Tot = NTA,
                  TRT.B.Tot = NTB,
                  Sum.Dis.Cured = SDC,
                  Cured = c,
                  Catch.Dis = New)
  
  # Predict S, I, R compartments by day
  while(day == -1 & d > 0 | (day != -1 & curr.day < day)){
    # if number of days is specified, calculate up until that day
    # else calculate until there is no more diseased individual
    
    # Update current day and c, u, d
    curr.day = curr.day + 1
    c = c + SDC
    u = d - SDC
    d = d + New - SDC
    
    # Calculate other variables
    h = N - d - c
    hs = round(g * h)
    s = hs + d
    a = s - u
    ns = round(PS * a)
    ad = round(PS * (d - u))
    as = ns - ad
    tps = rbinom(1, as, 1-tn)
    tpd = rbinom(1, ad, tp)
    tpt = tps + tpd
    NDAT = tpd + u
    NDTA = round(PTA * NDAT)
    NDTB = min(NDAT - NDTA, round(PTB * NDAT))
    DCA = min(NDTA, round(NDTA * A.eff))
    DCB = min(NDTB, round(NDTB * B.eff))
    if (curr.day == 2){
      NAT = tpt + tpt1
    }
    else{
      NAT = tpt + u
    }
    NTA = round(PTA * NAT)
    NTB = NAT - NTA
    SDC = DCA + DCB
    New = min(N - c - d, round(r * hs * d))
    
    data = c(curr.day, h, hs, d, s, a, u,
             ns, as, ad, tps, tpd, tpt,
             NDAT, NDTA, NDTA, DCA, DCB, NAT, NTA, NTB,
             SDC, c, New)
    df = rbind(df, data)
  }
  df$Cost = df$TRT.A.Tot * A.cost + df$TRT.B.Tot * B.cost + df$Screen.Tot * Test.cost
  return(df)
}

model.plot = function(model){
  cost = sum(model$Cost)
  sick = sum(model$Diseased)
  peak = max(model$Diseased)
  ggplot(data  = model, aes(x = Day)) +
    geom_line(aes(y = Healthy, color = "Healthy")) +
    geom_line(aes(y = Diseased, color = "Diseased")) +
    geom_line(aes(y = Cured, color = "Cured")) +
    scale_color_discrete(name = "Compartment") +
    labs(y = "Population", title = paste("Complex SIR Model: infected = " , sick, ", peak = ", peak, ", cost = $", cost, sep = "")) +
    theme_minimal()
}

df = complex.model.rand()
df

model.plot(df)

Checking variability by running 100 replications of the model

line.plot = function(model){
  list(geom_point(data = model, aes(x = Day, y = Healthy, color = "Healthy"), size = 0.8),
       geom_point(data = model, aes(x = Day, y = Diseased, color = "Diseased"), size = 0.8),
       geom_point(data = model, aes(x = Day, y = Cured, color = "Cured"), size = 0.8),
       geom_line(data = model, aes(x = Day, y = Healthy, color = "Healthy"), size = 0.1),
       geom_line(data = model, aes(x = Day, y = Diseased, color = "Diseased"), size = 0.1),
       geom_line(data = model, aes(x = Day, y = Cured, color = "Cured"), size = 0.1))
}

model.rep.plot = function(reps){
  plot = ggplot() + labs(y = "Population", title = paste("Complex SIR Model with", reps, "reps")) + theme_minimal()
  for(i in 1:reps){
    df = complex.model.rand()
    plot = plot + line.plot(df)
  }
  plot = plot + scale_color_discrete(name = "Compartment")
  plot
  return(plot)
}

model.rep.plot(10)

model.rep.plot(100)

accumulate_by <- function(dat, var) {
  var <- lazyeval::f_eval(var, dat)
  lvls <- plotly:::getLevels(var)
  dats <- lapply(seq_along(lvls), function(x) {
    cbind(dat[var %in% lvls[seq(1, x)], ], frame = lvls[[x]])
  })
  dplyr::bind_rows(dats)
}
visual_data = select(df, c(Day, Healthy, Diseased, Cured)) %>% tidyr::gather(key = Compartment, value = Population, Healthy, Diseased, Cured) %>% accumulate_by(~Day)
plot_ly(data = visual_data, type = "scatter", mode = "lines",
        x = ~Day, y = ~Population, color = ~Compartment, frame = ~frame) %>% 
  layout(hovermode = 'compare') %>%
  animation_opts(frame = 100, transition = 0, redraw = FALSE, mode = "immediate") %>%
  animation_slider(hide = T) %>%
  animation_button(x = 1, xanchor = "right", y = -0.1, yanchor = "bottom", label = "View Plot")

`arrange_()` is deprecated as of dplyr 0.7.0.
Please use `arrange()` instead.
See vignette('programming') for more help
This warning is displayed once every 8 hours.
Call `lifecycle::last_warnings()` to see where this warning was generated.